Mon. Not. R. Astron. Soc. 0Q0,[Tlfl6l(2Q10) Printed 17 February 2010 (MN KTeX style file v2.2) 



On the diffusive propagation of warps in thin accretion discs 



o 
in 



o 



> 

cn 
a^ 

(N 

o 
o 



Giuseppe Lodato^'^*, Daniel J. Price^'^ 

^ Dipartimento di Fisica, Universita Degli Studi di Milano, Via Celoria, 16, Milano, 20133, Italy 
^ Isaac Newton Institute for Mathematical Studies, 20 Clarkson Road, Cambridge CB3 OEH 

^Centre for Stellar and Planetary Astrophysics, School of Mathematical Sciences, Monash University, Clayton 3800, Australia 
'^School of Physics, University of Exeter, Stacker Rd, Exeter EX4 4QL 



Submitted: Revised: Accepted: 



ABSTRACT 

In this paper we revisit the issue of the propagation of warps in thin and viscous accretion 
discs. In this regime warps are know to propagate diffusively, with a diffusion coefficient ap- 
proximately i nversely proportional to the disc viscosity. Previous numerical investigations of 
this problem (iLodato & Pringlell2007l) did not find a good agreement between the numerical 
results and the predictions of the analytic theories of warp propagation, both in the linear and 
in the non-linear case. Here, we take advantage of a new, low-memory and highly efficient 
Smoothed Particle Hydrodynamics (SPH) code to run a large set of very high resolution sim- 
ulations (up to 20 million SPH particles) of warp propagation, implementing an isotropic disc 
viscosity in different ways, to investigate the origin of the discrepancy between the theory 
and the numerical r esults. We identify the cause of the discrepancy in an incorrect calibration 
of disc viscosity in iLodato & Pringlg (,2007). Our new and improved analysis now shows a 
remarkable agreement with the analytic theory both in the linear and in the non-linear regime, 
in terms of warp diffusion coefficient and precession rate. It is worth noting that the resulting 
diffusion coefficient is inversely proportional to the disc viscosity only for small amplitude 
warps and small values of the disc a coefficient {a < 0.1). For non-linear warps, the diffusion 
coefficient is a function of both radius and time, and is significantly smaller than the stan- 
dard value. Warped accretion discs are present in many contexts, from protostellar discs to 
accretion discs around supermassive black holes. In all such cases, the exact value of the warp 
diffusion coefficient may strongly affect the evolution of the system and therefore its careful 
evaluation is critical in order to correctly estimate the system dynamics. 
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1 INTRODUCTION 

Warped accretion discs may occur across a wide variety of astro- 
physical systems, from the large scales of accretion discs around 
supermassive black holes (SMBH), down to the small scales of 
planet forming discs. 

Observationally, warps are found in galactic binary systems, 
such as the hyperaccreting X-ray bi nary SS433 (Ise gelman et al. 
|2006|) and the X-ray binary Her X-1 l lWiiers& Prhiglell 19991), and 
in se veral microquasar s, including GRO J1655-40 jMartin et al.l 
l2008h and V4641 Sgr jMartin et al.ll2008l) . On the much less en- 
ergetic side, a warped protostellar disc is found around the young 
star KH 15D (Chiang & Murray-Clay 2004). 

Warps are also found in the thin accretion discs in Ac- 
tive Galactic Nuclei (AGN), as in the ca se of NGC 4258 
jHerrnstein et al.l [19961 : IPapaloizou et al. I ll998l) . The dynamics of 
warped accretion can play a fundamental role in these cases, as it 
in turn regulates the spin history of the growing SMBH and, as a 



conseq uence, its very ability to grow rapidly l lKing & Pringlel2006l: 
iKing e t al. 2008). 

The torques which produce the warp can be very differ- 
ent. For protost ellar discs, they include tidal interactio ns with a 
companion star jLarwood et al.|[l99^ : iMartin et al.ll2009l) , and dy- 
namical effects during the formation of the disc, which might 
affect the r elative orientation of the stellar spin and the plane- 
tary orbits teate et al.l boid) . For accretion discs around black 
holes there are additional torques arising from the general rela- 
tivistic Lense-Thirrin g prec e ssion around a spinning black hole 
teardeen & Petterson 1975; Scheuer & Feiler 1996; Ki ng et al. 
[2005; Lodato & Prin gle 2006; Martin et al.. 2007; Perego et al. 
2009), and se lf-induced warping caused by radiation pressure 
l lPringi3[l99^ . Recently, some attention has also been given to 
the process of disc warping and black hol e spin alig nment in the 
case of supermassive black hole binaries (Dotti et al. 2009). In all 
such cases, the evolution of the system is strongly dependent on the 
speed at which warping disturbances can propagate in the disc. 

Analytic theories of warp propagation have been discussed 



2 Lodato & Price 



extensively in the pa st IPapaloizou & Pringle|[T983l; |Pringlelll99j : 
|Papaloizou&LiiJl9 95: Ogilvie 1999l l200(]h (see Section^. These 
theories predict that while for thick discs warps should propagate 
as dispersive waves, with a velocity of the order of half the sound 
speed in the disc, in the limit of thin and viscous discs the propaga- 
tion is diffusive, wit h a diffusion coefficient inver sely proportional 
to the disc viscosity jPapaloizou & Pringle|[T983h . Numerical sim- 
ulations of warp propaga tion in the thic k dis c case have been per- 
forme d by Nelson & Papaloizoul ( Il999h and lNelson & Papaloizoul 
( l2000h . 

Numerical simulations of warp propagation for thin and vis- 
cous discs are much more challenging, because in order to prop- 
erly catch the warp dynamics it is essential to accurately resolve 
the vertical structure of the disc, which for very thin discs can be 
difficult. A first attempt at testing the analytical theory with nu- 
merical simulations in the thin disc regime has been performed by 
[Lodato & Pringle (2007) (hereafter LP07), using SPH. The results 
of lLPO?! showed some unexpected results: while for large values of 
the disc viscosity the warp diffusion coefficient appeared to scale 
inversely with viscosity, as predicted analytically, such behaviour 
was not found at low viscosities, already for small warp ampli- 
tudes. In this case, the diffusion coefficient appeared to be much 
smaller t han t heoretically predicted, implying (as discussed exten- 
sively by iLPOTi) some additional dissipation. Additionally, the in- 
ternal precession induced by the warp and predicted analytically 
was found to be strongly dependent on the specific implementation 
of viscosity and was not found to match the theoretical expecta- 
tions^ 

IlpovI discuss different possible explanations for such dis- 
agreement. On the one hand, it is quite possible that the limited 
numerical resolution of their simulations might have affected their 
result s. On t he other hand, strong supersonic motions were found 
m the lLPOTi simulations, which might result into shocks in the re- 
sulting flow and thus provide the required additional dissipation. 

In this paper, we want to systematically address all the issues 
left open by IlPOTI by checking both the numerical aspects of the 
problem and the physical effects involved. 

With regards to the numerical aspects , first of all we have used 
a different SPH code with respect to LP07, therefore validating one 
code against the other. Secondly, we have checked numerical con- 
vergence by running simulations using 20 million particles, that is a 
factor of ten larger thanlLPOVl (note that such simulations are among 
the largest SPH simulations of accretion discs performed to date). 
Since some of the effects reported by LP07 appeared to depend on 
the viscosity formulation, we have here tested two different pos- 
sible implementations of disc viscosity. Finally, we have modified 
our analysis procedure, so as to obtain a more quantitative evalu- 
ation of the uncertainties in the measured paramete rs. In o rder to 
test the physical effects which might determine the IlPOtI resul ts, 
we have paid attention to shocks, which were argued by LPOTI to 
be responsible for the additional dissipation. We have checked the 
importance of shocks by running simulations with different levels 
of bulk viscosity — varied independently of the shear viscosity — 
which is directly connected with shock dissipation. 

The paper is organised as follows. In section|2]we discuss the 
basic features of the analytic theory of warp propagation in both 
the linear and non-linear regime. In section[3]we detail the numeri- 
cal method that we have used to simulate the system, including the 
different implementations of disc viscosity that we use. In section|4] 
we describe the procedure we have used to analyse our results and 
extract from the simulations the warp diffusion parameters. In sec- 
tion|5]we present and discuss our main results for the warp diffusion 



and precession in both the linear and non-linear regime. Finally, in 
section[6]we draw our conclusions. 



2 ANALYTIC THEORIES OF WARP PROPAGATION 



We consider here, as in lLPOTl the propagation of warps in thin Ke- 
plerian accretion discs, rotating with angular velocity 0(7?), with 
surface density S(_R) and angular momentum per unit area L(7?). 
Here R should be interpreted as a 'spherical' coordinate. The local 
direction of L can be oriented arbitrarily in space, and the unit vec- 
tor l(i?) = L(_R) / L[R) defines its direction. If the disc is rotating 
around a central point mass M, then its rotation is Keplerian, with 
n = ^JGM/Ri andL(i?) = E{R)VGMR. 

The disc is warped whenever the direction identified by 1 
changes with radius. The warp amplitude may be characterised us- 
ing the dimensionless parameter -i/;, where 



dl{R) 



dR 



(1) 



The disc thickness is H — Cs/Q,, where Cs is the sound speed, 
and is the scale over which density and pressure change in the local 
2 direction. The disc aspect ratio is H / R, and we shall assume that 
H/R < 1. 

We use the standard I Shakura & SunvaevI jl973h prescription 
for the disc viscosity v, assumed here to be a standard, isotropic, 
Navier-Stokes viscosity: 



ctCsH — aQ,H 



(2) 



Warping disturbances can propagate in accretion discs in two 
different regimes, depending on the relative importance of pres- 
sure forces and viscous forces. If the disc is sufficiently thick, 
such that H/R > a, then t he warp propagates as a dispersive 
wave jPapaloizou &Linl 19951) . The equations of motion for a wave 
in the case where the disc is Keplerian a nd nearly inviscid are 
dLubow & Ogilviell2000l : lLubow et alj2002h 



dG 

dR 



(3) 



(4) 



where G is the disc internal torque in the horizontal plane (only). 
Note that these equations are valid only in the linear approximation 
for small warps, and that no general non-linear theory for the wave- 
like regime exists as yet (but see Ogilvie 2006). 

Here, we are mostly interested in the case where the disc is 
thin and viscous, such that H/R < a. In this case, the warp propa- 
gates diffusively jPapaloizou & P ringle 1983j, and can be approxi- 
mately described by the equation dPringldl 19921) 



dL 

dt 



3__a_ 

RdR 
RdR 



R^d_ 
E dR 



>iEi?'''2)L 



1^2R^ 



dR 



1^1 

2 



(5) 



I d (I p,^, 91 



In this equation, the terms proportional to describe the standard 
viscous evolution of a thin and flat disc. For small amplitude warps, 
ui = V (and thus qi = a), but for large amplitudes v\ can be 
affected by the warp. The terms proportional to V2 arise whenever 
the disc is warped and \d\/dR\ 7^ 0. According to Equation lO the 
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warp diffuses with a diffusion coefficient 1^2 . By analogy with the 
viscosity prescription (Eq. ^) we can define a second parameter 
Q2 so that 



z/2 = a2CsH = a20.H 



(6) 



It is clear that the nature of the evolution of a warped accretion 
disc is determined mainly by the relative values of a and 02- In the 
case of small warp a mplitude, ij) <^ H/R, and fo r viscosity such 
that H/R < a < l. lPapaloizou & Pringld 119831) have found the 
following relation between the two coefficient a and 02'- 



Q2 



1 

2^' 



(7) 



and therefore that the warp diffusi on coefficient is inversely propor- 
tional to the size of the viscosity. lOgilvid ( Il999h (hereafter lo9^ 
extends these approximate analytic results by use of an asymptotic 
expansion in terms of the small quantity H/R, but retaining the as- 
sumption of an isotropic (Navier-Stokes) viscosity. By this means 
he is able to take account of larger values of a and In the limit 
of a small amplitude warp (ip <C 1). , 099. finds 



1^2 _ 1 4(1 + 7a2) 
1^1 



2a2 4 + ^2 ' 

which includes higher order corrections in a. also computes 
the relation between a and Q2 for an arbitrarily large warp ampli- 
tude, for which there is no simple analytical expression, but that 
can be computed numerically. In order to analyse our results we 
will also make use of these numerical relations (see Sec.|5}. 

Finally, it should be noted that the full non-linear theory of 

I099I also includes some precessiona l torques, whi ch are not ac- 
counted for in the diffusion model bv lPringQ ( Il992h . because they 
arise at higher order in a. In interpreting our results, in some cases, 
we have added such terms in our simple diffusion model, by adding 
a term on the right-hand side of equation (|5}, in the form of (see 
lOgilvidl 19991 ) 



dL 

dt 



(9) 



where we have introduced a third coefficient related to preces- 
sional effects (with a corresponding ag — vj /Q.H'^). The non- 
linear theory of ,099 also provides an expression for the depen- 
dence of Q 3 on a and on t/). In the limit of a small amplitude warp 
(ip <^ 1) and to l owes t non-zero order in a, the precession coeffi- 
cient is given by ( |099|P] 



as = 3/8. 



(10) 



Taking account of large values of a, but still for small amplitude 
warps, the expected relation is given by (cf. lOgilvTe & Dubu jzodll 
Eq. 12) 



Q3 = 



3(1 - 2q^) 
2(4 + q2) 



(11) 



For non-linear war p ampli tudes, higher order corrections in both a 
and ?/) are given bv l099l Pl 



^ Note that this is given incorrectly as 3/4 in lPapaloizou & Pringld jl983h 
and correspondingly in lLPOTl 

2 The reader should note that l099l uses a slightly different notation, defin- 



3 NUMERICAL METHOD 

We have performed a series of three-dimensional SPH simulations 
of warped discs similar, though more extensive than, those per- 
formed bv lLP07t SPH is a Lagrangian scheme for solving the equa- 
tions of hydrodynamics in which fluid quantities and their deriva- 
tives are computed o n a set of A'^ particl es that follow the fluid mo- 
tion (see |Pricel2004l or lMonaghanll2005l for recent reviews). 

In this pap er we have used the PH ANTOM code, developed 
by D. Price (see lPrice & FederTathll20 1(f for another recent appli- 
cation). PHANTOM is a low-memory, highly efficient SPH code 
optimised for studying non-self-gravitating problems. The code is 
made very efficient by using a simple neighbour finding scheme 
based on a fixed (in this case, cylindrical) grid and linked lists of 
particles. In particular, the absence of overheads associated with the 
tree-code for computing gravitational forces as well as other opti- 
mis ations means that th e code is sig nificantly more efficient than 
the iBenz et al. 



ployed bv lLP07 



T990h / lBatd ( ll995l) -derived code previously em- 



The initial aim of using PHANTOM was, given the concerns in 

ILPOTI to be able to employ a much higher resolution in the calcu- 
lations. Whilst in the end we found this unnecessary (see Sec. |5]l, 
we have instead used our increased computing abil ity to survey a 
much wider parameter space than that explored by ILPOVI includ- 
ing a wide range of viscosity parameters, two different viscosity 
formulations and three different warp amplitudes. 



3.1 Navier-Stokes equations 

In this paper we compute the evolution of a viscous accretion disc 
by solving the Navier-Stokes equations for a viscous, compressible 
hydrodynamic gas in an external potential, given by 



dp 
It 
dv/_ 
dt 



= -pV • V, 
1 dS'' 



pot 1 



(12) 
(13) 



where the potential corresponds to the gravitational force from a 
central star or black hole of mass M at the origin, i.e.. 



GM . 



and the stress tensor is given by the usual expression 



2 \ dv'= 



(14) 



(15) 



where rj and ( are the shear and bulk viscosity coefficients respec- 
tively. Note that the kinematic shear viscosity v is related to the 
shear viscosity by — rj/p and similarly one may define the 
volume viscosity ^„ = ("/p. In the case of constant viscosity coef- 
ficients, the equations can be simplified to the vector form 



Itt 



VP 

P 



p \ 6 / p 



(16) 



The pressure is related to the density via a locally isothermal 
equation of state 



P = cl{R)p, 



(17) 



where the sound speed Cs is a prescribed function of (spherical) 
radius R = ^/a^^+j/^ip^z, given by 



ing the coefficients using Qi 



-3ai/2, Q2 = 012/2 and Qa = 03- 



Cs(i?) = Cs,0-R 



(18) 
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where q is given in Section [331 and the normalisation Cs,o deter- 
mines the d isc thickness. 

In the Ishakiira & SunvaevI ( Il973h a-parametrisation for the 
disc viscosity, the kinematic viscosity coefficient is given by Eq. 
l|2j. We consider two different methods for implementing Navier- 
Stokes viscosity in SPH, firstly based on a modification to the usual 
artificial viscosity term (similar to that employed by LP07) and sec- 
ondly based on a direct evaluation of the derivatives in dlSt and l ll3t 
respectively (see Sec. |3.2.4l l. The procedure used in the former case 
is described in Sec. 13. 2.3] setting r) and C, to zero in Eq.[T5] For the 
latter case we simply specify v in l llSt from the nominally input 
value for a using 



v{R) = a 



(19) 



where Cs (R) is specified according to Eq. l llSb and we assume a 
Keplerian rotation profile Q, = ^JGM/R?. The corrections to Ke- 
plerian rotation due to the pressure gradient are of order (H/R)^, 
which for the thin discs considered in this paper are very small 
(~ 10-*). 

3.2 SPH 

3.2.1 Hydrodynamics 

PHANTOM implements the fu ll variable smoot h ing l ength 
SPH formulation develop ed by IPrice & MonaghaiJ ( l2004h and 
IPrice & Monagh an (2007), whereby the smoothing length, h, and 
density, p, are mutually dependent via the density sum (for particle 
a) 



(20) 



Pa = y^^mtWab{ha), 
b 

which is an exact solution to ilH . and the relation 
/ \ 

ha = /itac , (21) 

where ma is the particle mass and Wgb = Vt^(|ra — rt\,ha) is 
the SPH smooth ing kernel (see e.g. lMonaghan|[l992l : |Pricell2004l : 
lMonagharill2005l for details). This results in a resolution that adapts 
to the local particle number density. Equations l |2Q| l and | |2U are 
iterated s elf-consistently using a Ne wton-Raphson method as de- 
scribed in lPrice & MonaghanI 1 20071) . where in this paper we have 
used ftfac ~ 1.2, giving approximately 58 neighbours per particle 
in a smooth distribution. 

The equations of motion il3\ take the form 



can be shown (see below) that the artificial viscosity corresponds 
directly to a Navier-Stokes type term and can thus be used, with 
minor adjustment of the parameters, to directly represent physical 
viscous diffusion (in doing so one would obviously discard the re- 
maining terms in Eq.|23] i.e., setting — rj = 0). The disadvantage 
of doing so is that the resultant viscosity coefficient consists of both 
shear and bulk components of viscosity, whereas for a disc the vis- 
cosity parameterisation l|2j should consist of shear viscosity only. 

The remaining part of the SPH algorithm is the time inte- 
gration algorithm, for which we use a standard leapfrog scheme 
equivalent to the velocity Verlet method. For efficiency we assign 
individual timesteps, set in factors of 2 from a nominal maximum 
timestep, such that only a subset of the particles is moved on the 
shortest timestep. With individual timesteps many of the conser- 
vation properties of the leapfrog algorithm are only approximately 
satisfied, however the scheme is significantly more efficient. 



3.2.2 Artificial viscosity 

The artific i al visc osity formulation in PHANTOM follows that of 
iMonagharil ( Il997h . with the averaging in the density and signal ve- 
locity changed slightly in order to more efficiently calculate the 
terms in i22l . We use 



^aa^ paVsig,a\'Vab ■ iabl, V^j, ' fab < 
3 Vab ■ fab > 



(24) 



where Vat = Va ~ Vb and the viscosity is only applied for ap- 
proaching particles (v^b ■ f^b < 0, i.e., converging flows). The 
signal velocity for hydrodynamics is given by 



-P |Vab ■ Tab] 



(25) 



where in general /? =2. The /? term in the signal velocity 
provides a non-linear term that was originally introduced to pre- 
vent particle pen etration in high Mach number shocks (see e.g. 
lMonaghan|[T98^ . 

For shock capturing — where the aim is to provide as little 
dissipation as possible whilst re solving shock str uctures — PHAN- 
TOM implements the .Morris & MonaghanI jl997h switch to reduce 
dissipation away from shocks, in which the dissipation parameter 



Q is evolved according to a source and decay equation 



dt 



^+Sa, 



Ta — ha/icrCs) 



(26) 



where a — 0.1, S = max(0, — V-v), Omin 
one would enforce amax ~ 1-0. 



0.05 and in general 



dt ^ ^ 



rrib 



^apl 



ViWab{ha) + 



■^ViWabihb) 
"bPb 



(22) 



where is a d imensionless quantity relat ed to the smoothing length 
gradients (see lPrice & MonaghanI I2OO7I for details) and the stress 
tensor is given by 



AV 



+ ria 



Pa 

dvl dvj 
dxi dxl 



+ 



C,a - -77c 



(23) 



The term in ( 



I is the artificial viscosity (discussed be- 
low) which is introduced in SPH in order to capture shocks and (to 
a lesser extent) to prevent interpenetration of particles. However, it 



3.2.3 Disc viscosity using the artificial viscosity term 

It has been known f or quite some time (e.g. lArtvmowicz & Lubowl 
ll994lMurravll996h that the artificial viscosity terms in SPH can be 
understood straightforwardly as numerical representations of sec- 
ond derivatives of the velocity. This is because the standard proce- 
dure for evaluating second derivatives in SPH is to use an integral 
formulation based on only the first derivative of the SPH kernel. 
For example the Laplacian for a scalar quantity A is represented by 
( lBrookshaw|[l985l : |Pricell2004 lMonagharj|200^ 



V^Aa 



rUb , 



—{Aa-Ab) 



fgb ■ VaVKgb 

|j"gb| 



(27) 



which is more clearly expressed by writing the kernel gradient as 

VaWab = VabFab, giviug 
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Figure 1. 3D stmcture of the waiped accretion disc from a representative 
20 million particle calculation in the large warp amplitude case (A = 0.5). 



(28) 



For a vector quanti t y the correspon ding expressions are 



tor a vector quanti t y the correspon c 
tespafiol & Revengall2003l : lMonaghaiJ|200^ 

^ \rab\ 



V(V-A) = -y!I!±\{S^^+2)iAai-rab)rab-Aai 



(29) 

Fab 



\rab\ 

(30) 

where 5^ =n i.e., the number of spatial dimensions. 

The above expressions mean that it is possible to give clear 
interpretation to the artificial viscosity terms, since for example in 
three dimensions we have, from {19\ and (|30t. 



Efnb , . ~ N Fab 
{Aab ■ rabJl T 

, pb \rab\ 



:V(V-A) + -V^A 



(31) 



For non-constant coeffici en ts the expressio ns are similar 
dClearv & MonaghanI 1 19991 : iMonaghad |2005|) . but with an 
average of the coefficients on the particles, for example 



V • (acVA) = - V — + ^b){Aa -At)^, 

^ pb \Tab\ 



(32) 



where IClearv & Monagh"ai] ( Il999h give alternative averaging pro- 
cedures more appropriate when the coefficients are discontinuous. 

Thus, for the artificial viscosity terms presented above (Sec, 
I3.2.2| l, we have 

AV 1 



Si; 



Vsig\rabl, 
1 AV II 

-a Vsig\rab\- 



(33) 
(34) 



Note that the factor |rai,| in the ^MonaghanI l ll997h formulation of 
viscosity used in PHANTOM differ s slightly from the factor h that 
would result from using the older IMonag hai] l ll992h formulation. 
The difference is only slight because b y definition within the ker- 
nel radius \rab\/h ^ 2, but use of the iMonaghanI i 19971) version 
avoids the need to account for divergences in the denominator when 
\rab\ ^ 0. 

In order to use th e artificial viscosity to represent a 
IShakura &Sunvaevl ( ll973l) disc viscosity, we therefore require sev- 
eral minor changes from the formulation appropriate for shocks 
given in Sec. 13.2.21 These changes are: 



(i) Viscosity should be applied for both approaching and reced- 
ing particles, 

(ii) The term in the signal velocity should be dropped such 
that Wsig = Cs, 

(iii) q^^ shou ld be multiplied by a factor h/\rab\, similar to the 
iMonag hai? /l99 2j) artificial visco s ity sc heme, and 

(iv) The Mor ris & MonaghM] ( 1 19971) switch should not be used 
i.e., should be treated as a constant. 

With these conditions the resultant 'artificial viscosity for a disc' is 
given by 



AV 1 

qa 



paCs,a\'Vab ■ rab\ 



Tab 



givmg 



AV 



>AV 



1 AV , 

-a c.h 

1 AV 7 

-a Csh. 
6 



(35) 



(36) 



(37) 



This is essentially the approach adopt ed bvlLP07land several earlier 
SPH accretion disc calculations (e.g. lArtvmowicz & Lubowil 19941 : 
Murray 1996), giving, from Eq. 



1 AvW 

To" IT- 



(38) 



where (h) is the azimuthally averaged (or, for a warped disc, shell 
averaged) smoothing length. The additional complication when us- 
ing a spatially variable smoothing length, addressed by LP07, is 
that in order to obtain a disc evolution corresponding to a single, 
uniform value of a thoughout the disc, it is necessary to setup the 
disc with a surface density profile such that {h)/H ~ const. This 
is discussed in Sec. 13. 31 below. 

The disadvantage of using the artificial viscosity term to rep- 
resent physical viscosity is that one inevitably ends up with a large 
and unwanted coefficient of bulk viscosity (Eq. l37l (. For a disc sim- 
ulation this is not so disadvantageous since in general V ■ v is not 
large, so although the coefficient is large, the term to which it is 
applied is small. However given that at least some of the deviations 
from the analytic theory found in .LP07i could possibly be explained 
by excess dissipation, it is desirable to perform simulations that ei- 
ther have no explicit bulk component or where the bulk viscosity is 
carefully controlled as would be the case when applied to shocks in 
the case of the usual artificial viscosity (Sec. l3.2."3t . 



3.2.4 Navier-Stokes viscosity implemented via two first 
derivatives 

A straightforward alternative to using the artificial viscosity term to 
model physical viscosity is simply to evaluate the terms in Eq. (123b 
direc tly. This is essentially the method proposed by Fleb be et al.l 
Il994l) . In this paper we evaluate the gradient terms using the stan- 
dard variable-smoothing-length gradient operator, given by 



dVa 1 ^ t i i\ 



pa^a 



dWabjha) 

dxi 



(39) 



where the coefficients of viscosity are set as discussed in Sec. 13. II 
Using this method one can in principle use zero bulk viscosity, by 
setting the bulk coefficient to zero and turning off any artificial vis- 
cosity terms. The danger with doing so is that any shocks that are 
present will not be treated appropriately and also that there is noth- 
ing to prevent particle interpenetration in the SPH scheme. Thus 
any such calculations should be treated with the appropriate degree 
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of caution. A better approach, and our default when using this for- 
mulation, is to set the physical bulk coefficient to zero in ( I23t , but 
to apply a small and carefully controlled amount of artificial vis- 
cosity to correctly dissipate shocks and approaching particles using 
the switches described in Sec. 13.2.21 The resulting coefficient of 
bulk viscosity in this case is, however, much lower than would be 
applied when using the artificial viscosity to mimic a disc viscosity 
(see Sec. 13.2.31 above). 



3.2.5 Navier-Stokes viscosity using direct second derivatives 

A further alternative, not considered in this paper, would be to di- 
rectly evaluate the second derivative terms resulting from the gra- 
dient of the stress tensor as in Eq. l |16t , using the standard rep- 
resentation of second derivatives in SPH given by Eqs. i29\ and 
l l3Qt . Indeed this f orms the basis of th e 'diss ipative particle dy- 
namics' scheme of lEspanol & Revengal ( |2003|) . The terms in this 
case are obviously similar to the formulation using artificial vis- 
cosity discussed above, except that the shear and bulk components 
can be set separately. The disadvantage is that the total angular 
momentum is no longer conserved because the dissipation is not 
applied along the line of sight joining the particles (meaning that 
mbVa X dva/dt 7^ 0). How serious a limitation this presents in 
practice for disc simulations has not been clarified, though it would 
be worthy of further investigation. 



3.2.6 Azimuthat averaging of SPH results 

In order to compare the results of the 3D SPH simulations with 
the simple diffusion equation l|5) it is necessary to compute 
azimuthally-averaged disc quantities from the SPH simulations. 
We perform this averaging by dividing the simulation domain uni- 
formly in R from Rin to i?out in A'^ = 350 spherical shells, with 
radial width A = (Rout — Rin)/N. The disc surface density in 
shell i is then given by the total mass divided by the disc surface 
corresponding to each shell, i.e., 



n[{R, + A/2)2 - {Ri - A/2)2 



(40) 



where the sum is performed over all particles in the shell. The av- 
erage angular momentum is computed as: 



Ni 



(41) 



where Ni is the number of particles in shell i. Finally, the local 
direction of the angular momentum vector can be computed using 



|Ji| 



(42) 



Examples of the resulting one-dimensional disc profiles are shown 
in Figures[3]|7]|9l[Tol[T2]and[T4] The above procedure for comput- 
ing the disc surface den sity has also been implemented as a feature 
in SPLASH ( |Pricell2007h . 



M = 1 (in code units). The gas particles are removed from the 
calculation inside a radius 7? — 0.5 (in code units). We distribute 
the particles using a Monte Carlo placement method such that the 
disc has a prescribed initial surface density profile, as described 
below. The particles are distributed in z so as to attain a Gaussian 
density profile in the vertical direction, with thickness H = Cs/il. 
The random particle placement, whilst simple, means that some 
settling of the disc occurs during the first few dynamical times of 
the simulation. 

The orbit of each particle is tilted such that the components of 
the unit vector 1 are given by 





A 
1 

A 



R — Rq 
R2 -Ri 



for R<Ri 



for Ri <R<R2 (43) 



for R> R2 



0, 



(44) 
(45) 



where Ri — 3.5 and R2 = 6.5 in code units, and _Ro ~ {Ri + 
R2) /2 = 5. The warp amplitude ■0 is then 



iP = R 



dR 



Rdh 
h dR 



(46) 



the maximum of which is attained at i? ~ i?o and is given by 



■kRo 



A 



2[Ri-R2) ^1- (^/2)2 



2.62- 



A 



(47) 



A three dimensional rendering of the resulting warped disc 
from one of our high resolution (20 million particle) calculations is 
shown in Figure[T] for the case of a high amplitude warp { A = 0.5). 
A cross section of the disc in a 2 million particle simulation with 
a low amplitude warp {A = 0.01) is shown in Figure |2l showing 
the slight bend induced in the dis c profi le. The initial shape of the 
warp is also plotted in Fig. 1 of lLPOTi and is shown by the lines 
corresponding to t = in Figures|7]and|9]of this paper. 

The disc extends from _Rin = 0.5 to i?out = 10, with a surface 
density profile, E, given by 



E(i?) = EoJ?" 




1 - \/ ^ 



(48) 



The parameter p is set to p = 3/2, as in IlPOTI such that, giving 
q — 3/4 in Eq. l llSt , the disc is uniformly resolved, in the sense 
that the smoothing length h cx P~^ ''^ is p roportional to the disc 
thickness H oc J?^'^*, as described in lLPOTt 

We choose the normalisation of the sound speed such that the 
aspect ratio of the disc at 7?o = 5 is H/R — 0.0133, corresponding 
to an aspect ratio at ii = 1 (in code units) of H/R = 0.02, in 
order to model a thin disc, where warps propagate primarily in the 
diffusive regime (see Section|2ll. 



3.3 Initial conditions 

Initial conditions are identical to those in IlpotI We use code units 
in which the gravitational constant G = 1, the central point mass 
M = 1 and the time unit is such that at a radius ii = 1 (in code 
units) the dynamical time is il^^ — 1. We place the gas particles in 
Keplerian orbits in the gravitational potential of a point mass with 



4 ANALYSIS 

4.1 ID Disc evolution 

We compare the time evolution of E and 1 from the SPH simulation 
with the one resulting from Equation (|5](. This is sol ved using stan- 
dar d finite difference methods, which are detailed in |Pringlelll992l) 
and IlpotI but with a different implementation of the zero torque 
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Figure 2. Cross section of the disc in the SPH calculations for a low amplitude warp (A = 0.01) at a resolution of 2 million particles. 



bound ary condition at the inner edge. While in |Pringlel l [l993) and 
IlpotI the zero torque condition is implemented in an approximate 
way, by artificially removing mass from the innermost cells in or- 
der to keep E close to zero, in this paper we directly enforce E = 
at the innermost cell. The two conditions are largely equivalent, 
but the shape of the surface density profile in the inner disc can be 
slightly modified. This is important because in turn it significantly 
affects the evaluation of a, as described in SectionlsTI 



4.2 Fitting procedure 

The main aim of this paper is to determine the relation between the 
two parameters a and 02, which describe the disc viscosity and the 
warp diffusion coefficient, respectively. In principle, the parameter 
a should be simply determined by the input viscosity coefficient 
in the SPH code as per Sec. 13.21 However, IlpotI did not find a 
perfect match between the nominal value of a as expected from the 
continuum limit of SPH and the a measured from the ID surface 
density evolution using Eq. ijSj, and therefore preferred to fit both 
param eters independently, to get the desired relation. In particular, 
IlpotI "stress that [they] do not perform an actual statistical fit of 
the viscosity coefficients, but simply choose them so as to match 
the evolution of the numerical simulation" . This point is discussed 
further in Section lSTI 

In this paper, we have implemented a statistical fitting proce- 
dure to check the calibration between the input a parameter and 
the value measured from comparing the SPH and the ID evolution 
of the disc using|5] The same procedure is further used to measure 
the warp diffusion coefficient 02 and the precession coefficient 03. 
This procedure is described below. Specific results for the calibra- 
tion of a and the estimate of 012 and 03 are reported in Sec. [5] 



4.2.1 Fittting for a 

The viscosity parameter a is primarily responsible for the evolution 
of the disc surface density E. In particular, given the shape of the 
surface density in our initial condition, the feature which is most 
directly related to a is the decline of the peak surface density in the 
inner disc. In order to obtain a we have therefore fitted the shape 
of the surface density profile close to the peak as resulting from 
the ID disc evolution at a given time to the SPH data. Thus, we 
have compared the data at t = 500 (in code units), and considered 
an annulus of radial width equal to 0.1 each side of the maximum. 
To obtain a, we have minimized the L2 norm of the difference 
between the ID evolution profile and the SPH data: 



= ^[E,-E^°(i?0]', 



(49) 



where the sum is taken over all shells (see previous section) within 
a radial distance 0.1 from the maximum. The value of T.™ at Ri 
is obtained by interpolation between the closest ID cells. The min- 
imum Ea is found by using a simple Newton-Raphson scheme. 
Starting from a trial value of a we iterate using: 



Ea{a„ + e) - Eaia„ - e) 



_Ea{a„ + e) + Ea{an - e) - 2Ea{an) 

where in the second line the first and the second derivatives of Ea 
are approximated by their finite difference value with respect to a 
small increment e. Once a minimum is found, we make sure that it 
is not a local minimum by checking that Ea is larger upon incre- 
menting a by 5e either side of the minimum. We also compute the 
1 — (J uncertainty on the fitted value of a by computing the distance 
from the minimum at which Ea is increased by a factor 2. 

An example of the best fit S profile compared to the averaged 
SPH profile is shown in Fig. [3] As it turns out, the best fit profile 
is in fact very close to the one computed using the input value for 
a, so this fitting procedure for a itself becomes unnecessary (see 
Sec. 15. II below). 

4.2.2 Fittting for 02 and as 

The values of Q2 and as are obtained using a similar procedure, 
though the computation of either requires, and is dependent on, an 
input value for the disc viscosity a — that is, using either the nom- 
inal or fitted a value. 

The diffusion coefficient Q2 is mostly responsible for the evo- 
lution of the profile of 1, and in particular, given our initial con- 
ditions, it affects the evolution of Ix around the warp radius at 
R — Ro. We therefore define the L2 norm of the difference 
between the ID evolution and the SPH evolution of Ix at time 
t — 1000 code units as: 



Ea2 = y^,[h,: 



(51) 



where the sum is taken over all shells within a radial distance equal 
to 3 from the warp radius. In practice, rather than fitting Q2 we fit 
the parameter /, defined as 02 = //(2q). The parameter / and its 
uncertainty are then obtained through minimization of Ea2 using a 
scheme analogous to the one used for a. 

The precession coefficient 03 affects primarily the evolution 
of ly around the warp radius. Its best fit value is thus obtained from 
the minimization of 



(52) 



(50) 



using a scheme analogous to the one used for a and 02, and de- 
pending on the input values of both of these. As for Q2 the sum 
is taken over all shells within a radial distance of 3 from the warp 
radius. 



5 RESULTS 

The initial aim of this paper was to perform simulations at a res- 
olution significantly higher than that employed bv lLPOTl in order 
to assess the effect of limited resolution. Having performed several 
calculations with 20 million SPH particles and finding results in- 
distinguishable from the lower resolution of 2 million particles, we 
have instead surveyed a wide range in parameter space, perform- 
ing a total of T8 simulations using 2 million particles together with 
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Figure 3. Left: surface density profile from tlie SPH simulation (solid black line) and from the ID evolution at t = and t = 500 in code units, for the case 
a = 0.3, using the input value (from Eg. 1381 (dashed green line) and the measured best fit value (dashed red line) of the disc viscosity for the calculation using 
artificial viscosity to model disc viscosity at an SPH resolution of 2 million particles. We find very close agreement between the input and fitted values. Right: 
Same, but using Navier-Stokes viscosity. In this case, the fit with the S profile is slightly better, giving smaller error bars in the fitted value of a. 



Series 


A 


Vise, 
type 


Bulk 
vise? 


Npart 


Symbols 


1 


0.01 


AY 


yes 


2 X 10^ 


red triangles 


2 


0.05 


AV 


yes 


2 X 10^ 


orange triangles 


3 


0.01 


AV 


yes 


2 X 10^ 


green triangles 


4 


0.05 


AV 


yes 


2 X 10^ 


cyan triangles 


5 


0.01 


NS 


no 


2 X lO** 


black squares 


6 


0.01 


NS 


switch 


2 X 10^ 


blue squares 


7 


0.5 


AV 


yes 


2 X 10^ 


magenta triangles 


8 


0.5 


AV 


yes 


2 X 10^ 


yellow triangles 



Table 1. Parameter settings for each of the 8 series of calculations per- 
formed in this paper, where each series consists of a set of simulations 
covering a range of disc viscosities a. The second column gives the ini- 
tial warp amplitude A used in Eq. j43K the third column shows whether the 
disc viscosity was represented using the modified artificial viscosity term 
(AV) or via a direct implementation of Navier-Stokes viscosity (NS). The 
fourth column shows whether or not bulk viscosity was applied (always true 
for the AV calculations). The resolution of the calculations is given in the 
fifth column and the symbols used to represent each series in Figures His] 
[6][8]and[T3]are given in the last column. 



the original 8 at 20 million particles. These consist of 8 series of 
simulations, each for a range of viscosity values. The parameters 
for each series are given in Table [T] where we have considered the 
effect of resolution (2 vs. 20 million particles), three different warp 
amplitudes {A = 0.01, 0.05 and 0.5), disc viscosity formulated 
either using the modified artificial viscosity (AV) or by a direct im- 
plementation of Navier Stokes terms (NS), considering the latter 
with zero bulk viscosity and subsequently with a small amount ap- 
plied using a switch. 



5.1 Calibration of the disc viscosity coefficient 

The best matching values of a (here referred to as Qflt) reported 
in Table 1 of ,LP07, do not follow the expected relation given by 
Eq. J38t . Although, for a given {h)/H, the relation between Qfit 
and is approximately linear, the slope is shallower than ex- 
pected. An accurate calibration of a is important, because in turn 
it is used as an input for the evaluation and fitting of the warp dif- 
fusion coefficient ct2. The disagreement found in IlpotI prompted 
us to examine the method used to calibrate a in greater detail, re- 
sulting in our implementation of the fitting procedure described in 
Section l42l— esse ntially a quantitative version of the procedure 
performed in lLPOTl 

In considering this issue, we have also explored the effect of 
the inner boundary condition of the ID disc evolution on the mea- 
surement of a. Indeed, the main feature which is used for the evalu- 
ation of a is the turnover of the surface density at small radii, which 
might well be affected by the specific boundary condition used. Us- 
ing the same condition employed by'LPO"^ (described in Sec. 14. II 
above), we found a similar relationship between Qflt and — 
that is, not matching the expectations from Eq. ( 1381 ). Furthermore, 
we found that the calculations using a Navier-Stokes viscosity also 
did not agree with the measured values. However, when the zero 
torque boundary condition was enforced exactly (see Sec. 14. It . the 
disagreement was removed. This is demonstrated in Fig. |4l which 
shows the best fit value Qflt versus the input value of a calculated 
using Eq. ([38}. Results are shown with triangles for the calcula- 
tions where the disc viscosity has been simulated using the SPH 
artificial viscosity, whilst squares correspond to calculations where 
Navier-Stokes viscosity terms have been implemented directly. All 
calculations are performed at a resolution of 2 million SPH par- 
ticles, except for the green and cyan triangles that use 20 million 
particles. The error bars show the \o errors from the fitting proce- 
dure. 
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Figure 4. Comparison between the input value of a and tlie measured value 
obtained fitting tlie SPH data to the ID evolution of E, with the zero torque 
boundary condition enforced exactly. Results are shown with triangles for 
the calculations where the disc viscosity has been simulated using the SPH 
artificial viscosity, whilst squares correspond to calculations where phys- 
ical viscosity terms have been implemented directly. All calculations are 
performed at a resolution of 2 million SPH particles, except for the green, 
cyan and yellow triangles that use 20 million particles. Error bars refer to 
the 1(7 errors from the fitting procedure used to obtain aat. 



ment was found at low a. Th e origin of the discrepancy between 
our results and those of lLPOTi Ues in the more accurate evaluation 
of the disc viscosity a discussed above and the resultant effect on 
the fit for Q2- 

The natural question is therefore: why do the new, more accu- 
rate results not agree with the standard theory? 

Initially, we will discuss only the simulations with a small 
warp amplitude A — 0.01. The effect of finite warp amplitude is 
discussed later in Section [5.2.5l 

5.2.1 Effect of resolution 

Our first attempt to reconcile the simulation with the theory was 
to check for resolution effects. In particular, these are thin discs, 
and the vertical scale-height is only moderately resolved (by only 
1.6 smoothing lengths) at a resolution of 2 million particles (see 
Fig. |2j. LP07 also mention the possibility of resolution effects in 
their simulations (which have the same resolution of 2 million par- 
ticles) but nevertheless demonstrate t hat the vertical density profile 
is very well reproduced (see Fig. 2 of IlpotI) . and hence argue that 
such effects should be negligible. We have therefore performed 6 
calculations at the higher resolution of 20 million particles (four 
with A — 0.01 and two with A = 0.05), for which the vertical 
scale-height is resolved with ~ 3.5 smoothing lengths. The results 
of these higher resolution simulations are shown with green and 
cyan triangles in Fig.[5]and show no significant difference with re- 
spect to the lower resolution case. Note that the calibration of a 
for these calculations, shown in Fig.|4]is very similar to the 2 mil- 
lion particle case, despite the very different smoothing lengths. We 
therefore conclude that resolution effects are not to blame for the 
discrepancy. 



As one can see from Fig. |4] the general agreement is very 
good. In other words, the procedure used by LP07 to simultane- 
ously fit a and a2 is no longer necessary and henceforth we simply 
fit the single parameter Q2 using the input value for a. 

Also worth noting from Fig. |4] is that, whilst the error bars 
ar smaller using the Navier-Stokes viscosity (squares) — due to 
an improved agreement of the profile of E near the inner boundary 
(see right panel of Fig. [3) — , there appears to be a small excess dis- 
sipation that occurs for low input a. In the case of Series 5 (black 
squares), this may be naturally attributed to the small amount of 
artificial viscosity we have applied. However, for Series 6 (blue 
squares), using zero artificial viscosity paradoxically results in an 
even greater dissipation. We attribute this to the increased random- 
ness of the particle distribution arising when no bulk viscosity is 
prese nt. A similar effect is seen in several other SPH simulations, 
e.g. in lPrice & Federratl] ( |20IO|) when /?^^ is too low. 

5.2 Results for the warp diffusion coefficient 

The results of the fitting procedure for the warp diffusion coeffi- 
cient 02 are shown in Fig. [5] for Series 1 to 5 discussed above. The 
solid line in the figure shows the simple relation 02 = l/(2a) ex- 
pected in the linear regime of warp propagation discussed in Sec- 
tion |2]([p^^^zou&^Wngl^[l98J). What is most striking is that 
the numerical results do not appear to match the 1/(2q) relation in 
almost any regime. In particular, at high a, the measured values of 
Q2 he much above the simple l/(2a) relation, while at low a they 
lie slightly below it. This is contrary to the result shown in iLPOl 
where agreement was found at high a, and a much larger disagree- 



5.2.2 Effect of viscosity formulation 

The results of lLP07 showed strong disagreement with the standard 
theory at low a. argue that a discrepancy between their re- 

sults and the standard theory at low a might arise because of en- 
hanced dissipation due to the presence of supersonic motions and 
hence shocks. With the recalibration of a discussed above the dis- 
agreement at low a is significantly reduced and, as we show later, 
can now be explained by the transition to the wave-like propagation 
regime. 

Nevertheless, when using the SPH artificial viscosity to rep- 
resent disc viscosity, one inevitably ends up with a large bulk vis- 
cosity coefficient (see Section |3.2.3| Eq. |37[l, wh ich could perhaps 
explain the excess dissipation invoked bv lLP07l In order to check 
whether bulk viscosity affects the magnitude of the warp diffusion 
coefficient, we have implemented a full Navier-Stokes viscosity, as 
described in Section[3] 

The results of the calculations using the Navier-Stokes viscos- 
ity and zero bulk viscosity and the artificial viscosity set to zero are 
shown in Fig.|6](blue squares), and indeed appear to show a signif- 
icant effect, lowering the measured values of Q2 at low a. But, as 
discussed in great detail in ILP071 this implies an even greater ex- 
cess dissipation. Thus, paradoxically, by removing bulk viscosity 
we have apparently increased the amount of dissipation present in 
the simulation. However, some caution is required here, since — as 
discussed in section [J.2.4| — one should be very careful about com- 
pletely removing all bulk viscosity from SPH simulations, since 
it plays a necessary role in preventing particle interpenetration, 
as well as providing physical dissipation in shocks, if shocks are 
present. 
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Figure 5. Relation between the warp diffusion coefficient 02 and tlie disc viscosity a, for warp amplitudes A = 0.01 (black squares, red and green triangles) 
and A = 0.05 (orange and cyan triangles). Squares use a direct implementation of the Navier-Stokes terms, whilst triangles use the SPH artificial viscosity to 
represent the disc viscosity. All calculati ons employ 2 million SPH p articles, except the cyan and green triangles, which use 20 million. The solid line shows 
the simple 0.2 = l/(2a) relation from lPapaloizou & Pringlj )l983l) . The long-dashed line includes the non-linear corrections due to finite values of a for 
small amplitude warps (Eq.[8). 



We have therefore also computed a series of simulations with 
Navier-Stokes shear viscosity, but with a small amount of artificial 
viscosity, applied only for approaching particles and using a switch 
to reduce its amplitude away from shocks (black squares in Figs.|5] 
and|6]l. Note that the effective bulk viscosity coefficient in this case 
is much smaller (at least a factor of 5) that that present in the sim- 
ulations that use artificial viscosity to represent disc viscosity (i.e. 
the red triangles in Fig.|5]l- Despite the very small change in the vis- 
cosity formulation, the results of this series of simulations (series 5, 
black squares in Fig.[5]l are in very good agreement with those pre- 
sented previously (series 1, red triangles). This leads us to conclude 
that the dramatic effect produced in Fig. [6] by removing all of the 
bulk viscosity, is most likely a numerical artefact. This is strength- 
ened by the results already discussed in Fig. |4l also showing that 
the simulations with zero bulk viscosity show a higher dissipation, 
that can be attributed to the increased randomness of the particle 
distribution when no bulk viscosity is present. 



5.2.3 Are we looking at the right theory? 

Having investigated the effect of both resolution and viscosity for- 
mulation, we may dispense with the possibility that the disagree- 
ment between the SPH simulations and the simple linear theory is 



simply a numerical artefact. What is worse, the recalibr ation o f a 
discussed in Section lSTI puts us in a worse position than lLPOTl as 
we now disagree with the line ar theo ry at both low and high values 
of a, contrary to the results of lLPOTl where the results appeared to 
agree at high a. In our case, the results are also more significant, 
because we have much bet ter sta tistics (86 simulations compared 
to the 10 shown in Fig. 8 o spanning a wider range in all of 

the parameters a, 02 and the warp amplitude. 

Clearly, though, the a values we are using here at not small, 
and therefore it may well be expected that the nominal 1 / {2a) re- 
lation is not appropriate, particularly at the high-a end of our pa- 
rameter range. Indeed, when we plot the relation between 0.2 and 
a predicted for finite values of a by Eq. ([8} (long-dashed line in 
Figs.[5]and[6]l, we recover an excellent agreement with the SPH 
results for a > 0.15 for low-amplitude warps. 

We thus conclude that there is no disagreement between our 
numerical results at high a and the linear theory of warp propaga- 
tion, once the effects of finite a are appropriately accounted for in 
the theory. The only remaining issue is a small disagreement at low 
a for s mall am plitude warps — though much smaller than the one 
found in lLPOTl which we discuss in the next section. 
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Figure 6. As in Fig. [5] but comparing only tlie two series using Navier- 
Stokes viscosity, witli (black suqares) and without (blue squares) a small 
amount of artificial bulk viscosity. Despite the nominally lower input vis- 
cosity, the calculations with zero bulk viscosity paradoxically show a larger 
dissipation, as evidenced by the lower magnitude of the warp diffusion co- 
efficient a2 . We attribute this to the danger of using zero bulk viscosity in 
SPH. 



5.2.4 Low viscosity behaviour and transition to the wave-like 
propagation regime 

In principle, there could be two explanations for the remaining 
small disagreement that we find for low viscosity: a numerical one, 
related to our procedure to fit the value of the diffusion coefficient 
02, and a physical one, related to an actual transition to a different 
propagation regime, such as the wave-like propagation regime (sec. 

There are good reasons to believe that both effects might play 
a role for small values of the viscosity coefficient. From the numer- 
ical point of view, we note that our fitting procedure is based on 
matching the solution of the simple diffusion equation to the SPH 
results at a given time, t = 1000 in code units, which corresponds 
to 0.4q in units of the viscous time at J? = 1. The viscous time is 
not only a measurement of the time needed for the overall viscous 
evolution of the disc to take place, but also of the time needed to 
smooth out the discreteness of the particle distribution in the ini- 
tial condition. We therefore might expect that, for small a, the SPH 
simulation at t = 1000 code units is still somewhat noisy, and that 
it might affect the evaluation of 02- That this is the case can already 
be seen from the fact that the error bars on 02 get larger at small 
a. Additionally, we have also found that for small ct the results are 
somewhat sensitive to the time at which the fit is performed. If the 
fit is performed at earlier times, the resulting value of 02 tends to 
be systematically shifted down by a small amount. 

The results of the calculation for low a might also be af- 
fected by the fact that warp propagation undergoes a transition to 
the wave-like regime (sec.|2j, which, for small warp amplitudes, is 
expected to occur at around a < 2-kH/\ « 0.07, where A is the 
wavelength of the perturbation (which in our case is A --^ 6 in code 
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Figure 7. Evolution of the shell averaged profiles of Ix from the SPH sim- 
ulations at 4 = and t = 500 in code units (sohd black lines), for the 
low viscosity case a = 0.065 and A = 0.01, compared to the results of 
a ID evolution assuming wave-like propagation with dissipation at the ap- 
propriate value of a, from Eqs. (3) and (4) (dashed red lines). The excellent 
agreement explains the deviations at low a seen in Fig. |5] as been due to 
the transition to the wave-like propagation regime at low viscosity. 



units). In order to test this hypothesis we compare the evolution of 
the SPH simulation to the ID evolution appropriate to the wave- 
like regime (Eqs.[3]and|4j, including dissipation corresponding to 
the input value of ct from the simulation. The results of this test are 
shown in Fig. |7] for a low amplitude warp with a = 0.065, from 
Series 1. The plot shows the profile of the a;-component of the unit 
vector 1 from the SPH simulation at t = and t = 500 in code 
units (solid black lines) compared to the ID evolution from Eqs. 
l[3j and Q (dashed red lines). The fact that the profiles agree well 
indicates that wave-like propagation, albeit with diffusion, can ex- 
plain the deviation from the purely diffusive propagation we have 
assumed when comparing with Eq. Q or Eq. l[8j. 

Given that there is no longer any disagreement between theory 
and the results of the SPH simulations, contrary to LP0 7. it is no 
longer necessary to invoke any additional dissipation involved in 
warp propagation. 

5.2.5 Large amplitude warps 

Having gained some confidence that the linear theory of warp prop- 
agation explains satisfactorily the evolution of the SPH simulations 
at low warp amplitudes, we may turn our attention to non-linear 
effects. 

Initially we have considered a small increase in the initial warp 
amplitude \.o A = 0.05 in Eq. i43t . The results of these calculations 
(Series 2) are shown with orange triangles in Fig. [5] and are essen- 
tially indistinguishable from the A = 0.01 case (Series 1 and 5), 
demonstrating that the linear theory is applicable also in this case. 

The resulting values of 02 from our fitting procedure for the 
A — 0.5 case (Series 7) are shown in Fig. [8] with magenta trian- 
gles, compared to the corresponding Series 1 for A — 0.01 (red 
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Figure 8. Relation between the warp diffusion coeiScient Q2 and a for 
large amplitude warps. The magenta triangles show the results for A = 
0.5, while the red triangles, for comparison, indicate the small amplitude 
{A = 0.01) case. The solid and long-dashed lines show the 1 / (2a) relation 
and the corrections expected for finite a, respectively. The short-dashed 
line shows the expected relation based on the non-hnear theory of l099l 
assuming a fixed ip = 0.55 (roughly comparable to the average ip value 
using A = 0.5). The large error bars in the high amplitude case at low 
a are because the non-linear warp evolution is not well fitted by a single 
value for a2 — physically manifested as a steepening of the warp profile 
observed in the simulations (see Fig. 110) that is not captured by the fitted 
linear profiles. 



triangles). One can immediately see that the fitted values of ot2, for 
small Q, are much smaller than the low amplitude case, and are 
characterized by a much larger uncertainty. 

At very large warp amplitudes, the linear theory of warp prop- 
agation in the diffusive regime is not applicable, and it is there- 
fore not surprising that, indeed, the SPH simulations in this regime 
{A = 0.5) do not match either the 1/(2q) relation or the more 
complete relation of Eq. ^ whi ch is non-linear in a but assumes 
linearity in the warp amplitude. l099l presents a non-linear theory 
of warp propagation for warps of any amplitude. We are now in a 
position to test this theory numerically. 

The complicating factor is that, in the 099 theory, 02 is a func- 
tion of the warp amplitude tp (related to A by Eq.|46l with the max- 
imum value given by Eq. 147 1 . We thus cannot associate a single 
value of a2 to each large amplitude simulation, as i/; is a function 
of radius (Eq.l46b and furthermore decreases as a function of time 
as the warp diffuses. Indeed, this is the reason for the large error 
bars in Fig.[8]when attempting to fit a single 02 value to the Ix pro- 
file, particularly at low a where in the l099l theory the dependency 
of 02 on is much stronger. What is possible is to check whether 
the deviations expected from such non-linear theory follow the ob- 
served trend, i.e. a general decrease value of the diffusion coeffi- 
cient. The short-dashed line in Fig. [8| shows the relation between 
02 and a based on the theory of |093 for a fixed value of ^ = 0.55, 
computed numerically based on a routine kindly provided to us by 
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Figure 9. Shell averaged profiles of Ix from the SPH simulation (solid black 
lines) at t = and t = 1000 (in code units), compared to the coitc- 
sponding profiles from the diffusive evolution, for the high viscosity case 
a = 0.43 and a strongly non-linear warp amplitude (A = 0.5), assuming 
the w arp diffusion coefficient to be a function of a and i/i, as predicted by 
I099I (dashed red lines). 



Gordon Ogilvie. We can indeed see that the non-linear theory does 
reproduce qualitatively the observed trend. 

As an additional test of the theory, we have also evolved the 
standard equation for diffusive evolution (Eq. |5), where the diffu- 
sion coefficients a\ (now different from a) and 02 are computed 
at each radius and at each timestep, based on Ogilvie's non-linear 
theory. In this case, we do not have any free parameter to fit: a is 
the nominal value of the viscosity coefficient, and a\ and a2 are 
prescribed functions of a and 1/; which we comput e at ea ch radius 
and time by evaluation of the relevant integrals in |099| using the 
routine provided. The resulting profiles of are shown in Fig. |9] 
for the case where a — 0.43. The black solid lines show the results 
of the SPH simulations at t = and t = 1000 (in code units), while 
the red dashed line show the corresponding profiles computed from 
Eq. The agreement of the non-linear theory with the SPH re- 
sults i s exce llent. Note that the warp diffusion coefficient in the the- 
ory of'099' depends also on the amount of bulk viscosity present in 
the disc. Indeed, in order to get the good match shown in Fig. |9] 
we have computed the warp diffusion coefficient assuming a bulk 
viscosity with magnitude Qb = 5a/ 3, as expected from the SPH 
formalism (see sect. 13. 2. 2t . 

At low disc viscosities (small a), the simulations show a 
steepening effect in the warp profile — the physical result of dif- 
ferent parts of the warp propagating at different speeds. We find 
that the steepening becomes stronger as the disc viscosity is re- 
duced (evident from the increase in the error bars seen in Fig.iJas 
Q — >■ 0) and in the limit of extremely small viscosity can result in 
a near-complete break in the disc. This is demonstrated in Fig. 1101 
which shows the equivalent of Fig. |9] for a low viscosity calcula- 
tion (a = 0.03) shown at intervals of t — 500 in code units and 
evolved as far as t = 3500. A three dimensional rendering of the 
disc structure is shown in Fig.[TT]at t — 1500 in code units. De- 
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Figure 11. Resulting 3D disc structure from the simulation shown in Fig. 1 101 
with a large amplitude warp in a low viscosity disc (a = 0.03), shown at 
t = 1500 (in code units). The steepening of the warp profile in this case 
results in a nearly complete break in the disc. At the employed resolution 
of 20 million particles we are able to resolve the thin but steady stream 
of material that is nevertheless still being transported inwards across the 
discontinuity (just visible in the Figure, as indicated by the arrow). 



Figure 10. Shell averaged profiles of Ix from the SPH simulation (black 
lines) at intervals of At = 500 up to t = 3500 (in code units), for the 
low viscosity case a = 0.03, using 20 million particles and a strongly non- 
linear warp amplitude {A = 0.5). The large warp amplitude simulations at 
low a show a strong steepening effect in the warp profile because different 
parts of the warp propagate at different speeds. A comparison with the non- 
linear theory is more difficult in this case because the predicted diffusion 
parameter for the disc can become negative, causing unphysical oscillations 
in the one-dimensional code. 



spite the apparent break the calculations are able to resolve a thin 
but steady stream of material continuing to accrete across the dis- 
continuity, indicating that the disc regions remain connected even 
in this extreme case. 

The simulation results in this case are more difficult to com- 
pare with theory because of the development of large unphysical 
oscillations around the warp in the one dimensional code. This is 
not unexpected from the non-linear theory, since in this regime the 
evolution of the di sc is not well described by a diffusion equation. 
In particular. [099I (Eq. 141) shows that 

_L:^ 



ai 



(53) 



and thus for large -tp (more specifically, when tp > v24q, which is 
the case here) the effective viscosity becomes negative, suggesting 
that diffusion is no longer an appropriate description. Whilst for 
small amplitude warps one expects a transition to the wave-like 
regime at low a (based on the dispersion relation derived from the 
linear wave equations, Eqs. l[3)-l|4}), it is not clear from theory at 
which value of a the transition occurs for large amplitude warps, 
and indeed whether such a transition occurs at all if tp is sufficiently 
la rge (Ogilvie, pr ivate communication). Furthermore, calculations 
bv logilvid lliooa) based on the evolution of a one-dimensional non- 
linear Schrodinger equation for inviscid, Keplerian discs suggest 
that wave-like propagation should not have a steepening effect for 
an isothermal equation of state (more specifically, if the adiabatic 
index 7 < 3), suggesting that the steepening we observe is indeed a 
diffusive rather than wave-like effect, though it is hard to be certain 
in the absence of a full non-linear theory for wave-like propagation. 



5.3 Precession 

The last remaining piece of the theory which needs to be tested is 
the one related to precession. As mentioned in sect. |2l the full the- 
ory of warp propagation presented in l099l predicts also the presence 
of internal precessional torques. These can be accounted for within 
the simple diffusion model of Pringle ( 1992) (Eq.O by adding an 
appropriate additional term (Eq. |9). We have included such terms 
and fitted the value of the corresponding new parameter as to the 
SPH data. In this case, the disc property which relates to precession 
and that we use to perform the fit is the profile of the j/-component 
of the angular momentum in the disc ly. 

Figure [12] shows the profiles of ly at t = 1000 code units 
(solid black lines) and the corresponding solution of the diffusion 
equation with added precessional term, using the fitted value of Q3 
(dashed red lines). The left panel refers to a = 0.43 and an SPH 
resolution of 2 million particles, while the right plot refers to a = 
0.46 and an SPH resolution of 20 million particles. Note that the 
amplitude of the induced precession — and therefore of the value 
of ly — is small, and therefore the SPH data are quite noisy for 
ly, especially at low SPH resolution. This implies a rather large 
uncertainty in the fitted value of 03. 

The resulting values of Q3 as a function of a are plotted in 
Figure [T3] for Series 1 to 4 (that is, for the small amplitude cases 
with disc viscosity modelled through the artificial viscosity term). 
Given the large uncertainty at the lower resolution and lowest warp 
amplitude (A = 0.01), error bars are shown only for the A — 0.05 
and higher resolution calculations for clarity (by way of compar- 
ison the error bars for the low resolution A — 0.01 calculations 
are roughly a factor of two larger than for the corresponding low 
resolution A — 0.05 results). The solid line shows the small am- 
plitude and a <^ 1 approximation as = 3/8 (Eq.llOt, whilst the 
long-dashed line gives the expected relation for small warp ampli- 
tudes but to higher order in a (Eq.[n)- Once again, provided the 
corrections for finite a are accounted for, we obtain a very good 
agreement between our numerical results and the theory, except for 
the lowest values of a. The deviations seen at low a are most likely 
due to the effect of bulk viscosity in the code which affects the ly 
profile more strongly than either Ix or l^ (simply due to the low 
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Figure 12. Left: Profile of ly aXt = 1000 (in code units) for the a = 0.43 case. The solid black line shows the SPH results at a resolution of 2 million particles. 
The dashed red line shows the expected profile from the diffusion plus precession equation with the fitted value of a^. Right: Profile of ly at t = 1000 code 
units for the a = 0.46 case, at an SPH resolution of 20 million particles (lines are the same as for the left panel). 



amplitude of ly), and is more significant when the disc viscosity is 
low. 

We note briefly that — by contrast with our earlier results 
for both a and 02 — the calculations utilising the Navier-Stokes 
implementation of viscosity (Series 5 and 6) show a strong dis- 
agreement with both the modified artificial viscosity calculations 
and the theoretical curves - the precession even reversing direction 
for a < 0.07. The fits for series 6 — i.e., the calculations that 
show good agreement in the 02 fits — show the fit to 03 rise with 
a (from negative values) and then flatten to around 03 ~ 0.23 at 
a > 0.2, in contrast to the results shown in Fig. [13] The results 
for series 5 show a similar trend but with much lower fitted val- 
ues. The errors to the fits are also significantly larger than for the 
artificial viscosity calculations. Whilst we can only speculate as to 
the reason for this disagreement, most likely it is an indication that 
higher resolution is needed (compared to using the modified artifi- 
cial viscosity) to evaluate the nested first derivatives (Eqs. 1391 and 
|13l l to the ac curacy required in order to measure the precessional 
contribution. IlPOtI also found the precession in their simulations to 
depend strongly on details of the viscosity formulation. 

Finally, let us consider the internal precession for large ampli- 
tude warps (A = 0.5, Series 7). In this case, as for the evaluation of 
the warp diffusion coefficient, we cannot simply associate a single 
value of 03 to our simulation, as it depends on the instantaneous 
value of 7/;, which is a function of R and t. However, we can still 
compare the profile of ly at a given time to the one expected from 
the non-linear theory of |099| . Once again, we stress that in this 
comparison we have not fitted any parameters, as the value of a is 
simply the input value in the simulation, while both and 03 are 
a prescribed function (obtained from ,099) of a, Qb = 5a/ 3 and 
if). The profile of ly for a = 0.43 and A = 0.5 is shown in Fig. 1141 
at t = and t = 500 (in code units). The solid black lines show 
the results of the SPH simulations while the dashed red lines re- 
fer to the solution of the diffusion equation with added precession. 



where th e coe fficients are computed directly from the non linear 
theory of 1 0991 . Note that, while in this large amplitude case the re- 
sulting shape of ly is a more complicated function than a simple 
oscillating function (as in the small amplitude case), the profile is 
reproduced surprisingly well by the 099 theory. To emphasize the 
importance of non-linear effects in this case, we also show with the 
dotted black line in Fig. [14] the profile of ly at t = 500 obtained 
from the ID evolution code neglecting the effects of non-linearity 
and simply adopting a constant 03 — 0.22, that is the value of the 
precession coefficient for a — 0.43 in the small amplitude limit. 
One can thus clearly see that the non-linearity in the determination 
of Qf3 is essential in order to reproduce the correct precession of the 
disc. 



6 CONCLUSIONS 

In this paper we have numerically tested the non-linear propaga- 
tion of warps in thin and viscous accretion discs. To this end, we 
have run very high resolution SPH simula tions o f warped accre- 
tion discs, extending the previous work of lLP07l to cover a much 
wider region of the parameter space. In some simul ations , we have 
increased the numerical resolution with respect to ILP07| by using 
ten times as many particles. We have also checked the effect of two 
different implementation of the disc viscosity. 

O ur new and improved results correct upon the previous re- 
sults of lLP07l who had found a disagreement between their simu- 
lations and the analytical theories of warp propagation. On the con- 
trary, our results are in specta cular agreement with the non-linear 
theory of warp propagation of l099l Some specific features of this 
theory, confirmed by our simulations, are worth recalling: 

(i) For moderate values of q > 0.1, the warp diffusion coef- 
ficient V2 is not proportional to 1/a, where a is the disc viscos- 
ity coefficient, but follows the slightly more complex relation, Eq. 
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Figure 13. Relation between the precession coefficient 03 and the disc vis- 
cosity a, for warp amplitudes A = 0.01 (red and green triangles) and 
A = 0.05 (orange and cyan triangles). All calculations employ 2 million 
SPH particles, except the cyan and green triangles, which use 20 million. 
The solid line shows the expected precession rate in the limit of small a, 
while the dashed line shows the relation betw een 03 and a expected for 
small amplitude warps from the theory of l099l (Eg. lilt . 

((Sj. The 'standard' 1/a behaviour is only recovered for smaller 
values of a and for small amplitude warps. Note that a value of 
Q « 0.1 is expected ba sed on observations of accreting binary sys- 
tems ( iKing et alj|2007[) . 

(ii) For large amplitude warps, the relation between V2 and v is 
much flatter than the 1/a relation, corresponding to a more uni- 
form (with respect to the disc viscosity) but also much less efficient 
diffusion of the warp at low a compared to the linear case. Our 
simulations, which are characterised by a warp amplitude i/; « 1 
are reasonably well described by an almost constant 02 ~ 2.5. 

(iii) In general, for non-linear warps, the warp diffusion coeffi- 
cient is a function of the warp amplitude, which is itself a function 
of position and time. For a proper calculation of the warp evolu- 
tion in a simple ID diffusion code it is essential to include such 
dependence. We stress that this can and should be done in any warp 
diffusion code. 

(iv) The non-linear theory also predicts the appearance of inter- 
nal precessional torque s. Al so such torques are well described by 
the non-linear theory of l099l and can be easily included in ID mod- 
els by the simple addition of an extra term in the evolution equation, 
as discussed in the text. 

(v) For large warp amplitudes and small viscosity (1/; > \/24q) 
the evolution of the system is not well described by a simple diffu- 
sion equation and a full numerical approach is thus needed in such 
cases. 

All of the above aspects of warp propagation are reproduced 
faithfully in our simulations and can deeply affect the stmcture and 
evolution of the disc. We stress that their implementation (except 
for the last point) is relatively easy and recommended in simple ID 
models of warped accretion disc evolution. 

Before concluding, we should add a word of caution. All of 



Figure 14. Profile of ly axt = 500 time units for the case A = 0.5 and a = 
0.43. The solid black lines refer to the SPH simulations, while the dashed 
red line show the result of the evolution of the diffusion plus precession 
code, with non linear warp parameter 02 and 03 computed based on l099l 
theory of warp propagation. For comparison, we also show with the dotted 
black line the profile at t = 500 obtained from the simple ID evolution 
model assuming a constant 03 = 0.22, appropriate for this value of a in 
the linear regime (i/) <g; 1). 



the above results refer to warp propagation in a disc where vis- 
cosity is a standard Navier-Stokes viscosity, and in particular to 
the case where it is isotropic. Accretion disc viscosity is gener- 
ally thought to be due to turbulence driven b y some disc instability , 
such as the magneto-rotation al instabi lity (Balbus & Hawlev 199 1|) 
or gravitational instabilities jLodato & Ric e ^2004 ; Lodato 20Qi). 
In such cases, it is not obvious that the induced transport can be de- 
scribed in terms of an isotropic viscosity coefficient, and some of 
the abo ve results, in particular concerning the precessional terms 
( IlPOTI) . might be affected. 
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